/** Fichier implémentant une forme avec une surface définie par une équation du
 * 4ème degré.
 *
 * Ce fichier contient l'implémentation de la classe Tricurte.
 *
 * @file Tricurte.cpp
 * @author Johan "Solo" Luisier
 * @date 20 Février 2013
 */
#include "Tricurte.hpp"

#include "Algebre/Solveur.hpp"
#include "Utils/ATL.hpp"

using namespace std;

Tricurte::Tricurte()
    : Forme( "Tricurte" ), Alpha(), Beta(), Kappa()
{
}

Tricurte::~Tricurte()
{}

Coord Tricurte::trouveIntersection( const Rayon& ray ) const
{
    Coord pt( VersLocal( ray.origine() ) );
    Vecteur dir( VersLocal( ray.direction() ) );

  /**
   * Il faut résoudre une équation super moche du 4ème degré
   * @f$a\,t^4+b\,t^2+c\,t^2+d\,t+e=0@f$. On définit les coefficients
   * ici.
   */

    double a( dir.norme() - Beta * ( 2. * dir.x() * dir.x() - dir.z() * dir.z() ) * ( 2. * dir.y() * dir.y() - dir.z() * dir.z() ) );
    double b( 4. * Beta * ( Kappa * ( dir.x() * dir.x() - dir.y() * dir.y() ) * dir.y() -
                            dir.x() * ( 2. * dir.x() * pt.x() - dir.z() * pt.z() ) * pt.x() +
                            ( dir.y() * dir.y() * pt.z() + dir.y() * dir.z() * pt.y() - dir.z() * dir.z() * pt.z() ) +
                            4. * dir * pt * dir.norme() ) );
    double c( 2. * Alpha * puissance< 2u >( Kappa ) * dir.norme() +
              2. * Beta * ( puissance< 2u >( Kappa ) * dir.norme() +
                            2 * Kappa * ( dir.x() * dir.x() * pt.z() +
                                          2. * dir.x() * dir.z() * pt.x() -
                                          dir.y() * ( dir.y() * pt.z() +
                                                      2. * dir.z() * pt.y() ) ) -
                            dir.x() * dir.x() * ( 2. * pt.x() * pt.y() - pt.z() * pt.z() ) -
                            4. * dir.x() * ( 2. * dir.y() * pt.y() - dir.z() * pt.z() ) * pt.x() -
                            dir.y() * dir.y() * ( 2. * pt.x() * pt.x() - pt.z() * pt.z() ) +
                            4. * dir.y() * dir.z() * pt.y() * pt.z() +
                            dir.z() * dir.z() * ( pt.x() * pt.x() + pt.y() * pt.y() - 3. * pt.z() * pt.z() ) +
                            2. * ( dir.x() * dir.x() * ( 2. * pt.x() * pt.x() + pt.r2() ) +
                                   4. * dir.x() * ( dir.y() * pt.y() + dir.z() * pt.z() ) * pt.x() +
                                   dir.y() * dir.y() * ( 2. * pt.y() + pt.r2() ) +
                                   4. * dir.y() * dir.z() * pt.y() * pt.z() +
                                   dir.z() * dir.z() * ( 2. * pt.z() * pt.z() + pt.r2() ) ) ) );
    double d( 4. * Alpha * puissance< 2u >( Kappa ) * dir * pt +
              4. * Beta * ( puissance< 2u >( Kappa ) * dir * pt + Kappa * ( 2. * dir.x() * pt.x() * pt.z() -
                                                                            2. * dir.y() * pt.y() * pt.z() +
                                                                            dir.z() * ( pt.x() * pt.x() - pt.y() * pt.y() ) ) -
                            dir.x() * pt.x() * ( 2. * pt.y() * pt.y() - pt.z() * pt.z() ) -
                            dir.y() * ( 2. * pt.x() * pt.x() - pt.z() * pt.z() ) * pt.x() +
                            dir.z() * ( pt.r2() - 2. * pt.z() * pt.z() ) +
                            4. * ( dir * pt ) * pt.r2() ) );
    double e( puissance< 2u >( Alpha * Kappa * Kappa ) + 2. * Alpha * Kappa * pt.r2() -
              Beta * ( Kappa * Kappa - 2. * Kappa * pt.z() - 2. * pt.x() * pt.x() + pt.z() * pt.z() ) *
              ( Kappa *  Kappa + 2. * Kappa * pt.z() - 2. * pt.y() * pt.y() + pt.z() * pt.z() ) +
              pt.r2() );

    MoteurResolution -> setParams( e, d, c, b, a );

    vector< double > solutions;

    MoteurResolution -> resoud( solutions );

    vector< double >::const_iterator It( solutions.begin() ),
            End( solutions.end() );

    Coord solution;

    for ( ; It != End; It++ )
    {
        if ( *It > -1.e-7 )
        {
            solution = pt;
            solution += *It * dir;
            return VersGlobal( solution );
        }
    }

    return NullePart;
}

/** Calcule la valeur du champ de normale au point donné.
 *
 * On utilise la dérivée de l'équation de la surface au point donné.
 *
 * @warning Ne vérifie pas si le point est bien à la surface de l'objet.
 *
 * @param[in] pt point de la surface pour lequel on désire connaître la normale.
 *
 * @return la normale à l'ellipsoïde par le point @c pt.
 */
VecteurU Tricurte::normale( const Coord& pt ) const
{
    Coord local( VersLocal( pt ) );

    return VecteurU( 4. * puissance< 3u >( local.x() ) - 4. * local.x() * ( ( 2 * Beta - 1 ) * puissance< 2u >( local.y() ) - ( Beta + 1 ) * puissance< 2u >( local.z() ) - 2 * Beta * Kappa * local.z() - ( Alpha + Beta ) * puissance< 2u >( Kappa ) ),
                     4. * puissance< 3u >( local.y() )  - 4. * local.y() * ( ( 2 * Beta - 1 ) * puissance< 2u >( local.x() ) - ( Beta + 1 ) * puissance< 2u >( local.z() ) + 2 * Beta  * Kappa * local.z() - ( Alpha + Beta ) * puissance< 2u >( Kappa ) ),
                     -16. * ( ( Beta - 1 ) * puissance< 3u >( local.z() ) - local.z() * ( ( Beta + 1 ) * puissance< 2u >( local.x() ) + ( Beta + 1 ) * puissance< 2u >( local.y() ) + ( Alpha + Beta ) * puissance< 2u >( Kappa ) ) - Beta * Kappa * ( puissance< 2u >( local.x() ) - puissance< 2u >( local.y() ) ) ) );
}

/**
 * @todo Implémenter cette méthode.
 */
bool Tricurte::estDans( const Coord& pt ) const
{
    return false;
}

bool Tricurte::estSur( const Coord& pt ) const
{
    Coord local( VersLocal( pt ) );

    return abs( puissance< 2u >( local.r2() - Alpha * puissance< 2u >( Kappa ) ) - Beta * ( puissance< 2u >( local.z() - Kappa ) - 2. * puissance< 2u >( local.x() ) ) * ( puissance< 2u >( local.z() + Kappa ) - 2. * puissance< 2u >( local.y() ) ) ) < 1.e-7;
}

ostream& Tricurte::affiche( ostream& os ) const
{
    return os;
}

Coord Tricurte::intrinseque( const Coord& pt ) const
{
    if ( ! estSur( pt ) )
    {
        return NullePart;
    }

    return NullePart;
}
